Phase transition(s) in finite density QCD 
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The Grand Canonical formalism is generally used in numerical simulations of finite 
density QCD since it allows free mobility in the chemical potential /i. We show that 
special care has to be used in extracting numerical results to avoid dramatic rounding 
effects and spurious transition signals. If we analyze data correctly, with reasonable 
statistics, no signal of first order phase transition is present and results using the Glasgow 
prescription are practically coincident with the ones obtained using the modulus of the 
fermionic determinant. 

1. INTRODUCTION 

Finite density QCD is a field where the lattice technology is still at an early stage. 
Even if considerable progress has been achieved in the investigations of equation of state 
of QCD at finite temperature and zero chemical potential using the lattice approach, the 
present situation of these studies at non zero density is not so satisfactory. 

As well known, the complex nature of the determinant of the Dirac operator at finite 
chemical potential makes it impossible to use standard simulation algorithms based on 
positive-definite probability distribution functions. 

Theoretical predictions in the lattice formulation are mainly limited to the infinite 
coupling limit and rely on some approximation (e.g. mean field). In this limit they predict 
a strong first order saturation transition at a value for the chemical potential significantly 
smaller than one third of the baryon mass (as one could naively expect considering the 
quarks confined inside hadrons but ignoring their binding energy) . The inclusion of some 
(3 dependence in analytical calculations indicates that the mean field critical density and 
|m# converge toward a common limit in the physically relevant scaling region ||. 

Until recently only two results were available in numerical simulations, both in the 
strong coupling limit. The first is based on a representation of the partition function as 
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a system of monomers, dimers and baryonic loops (MDP) ||. The results show a clear 
first order signal at a critical chemical potential fi c only slightly larger then the mean field 
prediction. 

The second is based on an expansion in powers of the fugacity variable z = e M : the 
Grand Canonical Partition Function formalism [||]. The main advantage of this formal- 
ism is the free mobility in the chemical potential once we have all the coefficients of the 
expansion. The non positivity of the integration measure can be overcome generating 
configurations accordingly to the Dirac determinant at zero chemical potential and defin- 
ing Z = (det A(/i)/ det A(fi = 0)) (also known as Glasgow algorithm). The non positivity 
of the Dirac operator reflects in a non positivity of the fugacity coefficients, at least for 
reasonable statistics. Nevertheless the complex coefficients can be used to compute the 
Lee- Yang zeros in the fugacity plane and extract some physical information. Using this 
technique in the (3 = case, evidence has been found for a non saturation first order 
transition at a chemical potential in good agreement with the one extracted from MDP 
simulations Another interesting result of the Glasgow algorithm was the appearance 
of an early onset for the number density at a value fi D going to zero in the chiral limit as in 
simulations where the modulus of the fermionic determinant is considered. In the latter 
case it can be explained as the effect of light baryonic states formed from ordinary quarks 
and "conjugate" quarks interacting among themselves through the complex conjugate of 
the Dirac operator [Q. In QCD such states are certainly not present and the meaning of 
this early onset has never been understood. 

The aim of this paper is to show that the GCPF has some numerical subtleties to be 
considered. Once rounding effects are under control the data show no signals of phase 
transition With reasonable statistics the only structure left is the early onset in the 
number density and no signals of first order phase transition survive: numerical results 
are practically coincident with the theory defined using the modulus of the determinant. 

In the second section we introduce the GCPF formalism and show how (wrong) numer- 
ical manipulations can simulate a first order phase transition in the theory with quarks 
and conjugate quarks. The third section is devoted to the study of numerical instabilities 
in the evaluation of coefficients of fugacity expansion. In the fourth the correct numerical 
analysis is reported together with conclusions. 

2. MODULUS OF THE DETERMINANT 

The contributions of the modulus of the Dirac determinant (A) and its phase </>a can 
be separated as 



where (e t ^ A )\\ accounts for the mean value of the cosine of the phase computed with the 
probability distribution function of the pure gauge theory times | det A|. It is clear that, 
in the thermodynamic limit, (e l ^ A )\\ gives a net contribution to the free energy density 
only in the case in which it falls off exponentially with the lattice volume. This is not what 
happens in + 1 dimensional QCD where the average of the cosine of the phase is equal 
to 1 0. To check if this is the case in 3 + 1 dimensions we performed simulations of QCD 
at (3 = on different lattices using the modulus of the Dirac determinant f| . We decided 
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to generate configurations randomly after checking that the effective fermionic action is a 
flat function of the pure gauge energy around E = 1. The GCPF was used to reconstruct 
the coefficients of the fugacity expansion from the eigenvalues of the propagator matrix 

3V 

det A(m q , fj) = z 3V det (P{m q ) - z^l) = ^ a n {m q )z n (2) 

n=-3V 

where V is the lattice volume. 

We considered the number density n(ju) as well as its derivative that, in the thermody- 
namic limit, is proportional to the radial distribution of Lee- Yang zeros in the complex 
fugacity plane p(|e^|): 

An( M ) = 47re 2 M|e1). (3) 

The plots show results in some way contrary to common wisdom: if the system spatial 
extent is large enough first order phase transition signals are evident (fig. [I]). The small \x 
behaviour is characterized by an onset threshold that goes to zero in the chiral limit (fig. 
H). For larger fi a clear structure develops for n(p, c ) ~ 0.6 as well as at the saturation 
point. At small and intermediate masses, contrary to MDP results, the first order critical 
point fi c does not bring the system directly to saturation; nevertheless its position is in 
good agreement with what found using the MDP algorithm. From \i = up to \i c there 
is also striking agreement with results obtained using the Glasgow algorithm. For /i > // c 
the agreement is only qualitative but both methods show a saturation transition at /i s . 
On the contrary, in the quenched case, no signals of transition have been found indicating 
that the unquenching is essential to reproduce the finite density physics [§|],f|. 

We have also used the Glasgow prescription (i.e. setting to zero the averaged a n with 
negative real part) and calculated the number density obtaining identical results. The 



agreement between our calculations and the Glasgow algorithm is surprising since the 
two models are believed to have completely different baryonic spectrum. Perhaps more 
surprising is the behaviour for large masses. In this limit we can expect, on phenomeno- 
logical grounds, a clear transition signal since all the nucleons become confined in a single 
spatial lattice site and do not interact each other. This is not what is observed and \i c 
disappears as soon as m q > 0.7 (similar behaviour is reported also in ||). 

To have an independent check for this result we decided to diagonalize directly the 
Dirac matrix A(m q = 0) for various values of /i: for large masses a clear transition signal 
appears in agreement with extrapolation of data from MDP (fig. [3]) ||. 

The discrepancy between results obtained from P and A eigenvalues points therefore 
to the existence of numerical problems in the evaluation of thermodynamical quantities. 

3. ROUNDING EFFECTS 

There are three possible sources of numerical problems: i) diagonalization of A, ii) 
diagonalization of P, Hi) determination of the fugacity expansion coefficients a n . 

We used a standard NAG library routine to perform the diagonalization of A and P. 
Using the eigenvalues of the two matrices (constructed from the same gauge configuration) 
we verified that the equality between det A and z 3V det(P — z~ l I) (the latter evaluated 
directly from P eigenvalues) holds in the whole range of p, and masses we considered. At 
p = 0, we also found a perfect agreement between the A eigenvalues computed with this 
routine with the ones obtained using a standard Lanczos algorithm. 

We conclude that the diagonalization procedure is stable for any value of m q and p 
and that the numerical problems can only be due to the manipulations necessary to go 
from the P eigenvalues to the GCPF coefficients a n . The only feasible way to perform 
this calculation is to adopt a standard recursion method. We use the first n eigenvalues 
to calculate the coefficients of a polynomial Q n of degree n: 

4 n) = afc? - X n at 1} k > 1 (4) 

where ajf ^ is the coefficient of order k of Q n -\ and {Aj} are the eigenvalues. 

This algorithm has problems similar to the deflating technique for finding the zeros of 
large polynomials ||10|| . We can get the right answer only if we are able to reproduce the 
(approximate) symmetries of the correct result at each intermediate step of calculation. In 
this case the coefficients never grow too much and we do not have to rely on cancellation 
of large numbers to get the right answer. 

To clarify this point it is useful to consider a simple case: a single configuration whose 
{Aj} are uniformly distributed on two circles of radius p and p _1 (to fulfill the A, 1/A* 
symmetry). The fugacity polynomial contains only three non vanishing terms: a±^v = 1 
and ao = p v + p~ v . In the infinite volume limit this corresponds to a first order saturation 
transition with the number density that can be approximated by a Heaviside 9 function. 

If we use (|p to calculate the coefficients we can get different results depending on the 
ordering of the eigenvalues ||. If the eigenvalues are phase-ordered, at any intermediate 
step we calculate the coefficients of a polynomial whose zeros lay on two arcs of circle 
of increasing length. These coefficients are non zero and of order 0(p n ): in this case 
rounding propagation easily prevents to obtain the correct answer. This happens already 




Figure 3. Number density in a 6 3 x 4 lat- 
tice at (3 — and m q = 1.5 with shuf- 
fled (solid line) and unshufHed (diamonds) 
eigenvalues; the same quantity for a 4 4 lat- 
tice from the A eigenvalues (squares). 



Figure 4. Number density in a 6 3 x 4 lat- 
tice at (5 — and m q = 0.1 using shuf- 
fled (solid line) and unshufHed (diamonds) 
eigenvalues. 



for relatively small V and forbids the symmetries to be realized in the final results. If 
we use randomly ordered eigenvalues we modify this scenario since the symmetries are 
(almost) enforced at each intermediate step as well as in the final result ||: the coefficients 
of Q n never grow too much and rounding effects are better under control. 

The eigenvalues from the diagonalization routine can be (partly) phase ordered (and 
this is what often happens with our NAG routine): they can not be considered a good 
input for the routine that evaluates the coefficients. Rounding effects become more and 
more probable as the volume increases since in this case the coefficients are bigger and a 
smaller fraction of ordered {Aj} can suffice to create problems. 

For actual configurations the eigenvalues are distributed in two (almost) circular strips 
in the complex plane. In this case we have to pay attention to the ordering of their 
modulus too. The shuffling procedure solve both problems since at each intermediate 
step the zeros of P n have, in the complex plane, the same distribution as the {Aj}. When 
we do it we get from the coefficients results indistinguishable from the ones obtained 
computing the fermionic determinant directly from the propagator matrix eigenvalues. In 
fig. ^ is shown that once we solve the numerical problems we can reconcile results from 
P and A eigenvalues. 



4. STRONG COUPLING RESULTS 

After we have removed the numerical artifacts we can reconsider the strong coupling 
results at small quark masses. The plot in fig. |] shows that data, if correctly analyzed, give 
no evidence of first order transition (the plot refers to data obtained using the Glasgow 
prescription but, using the modulus, we get the same results f§). No signal is present in 



the position indicated from MDP results and the only structure left is the onset fi , once 
again consistent with the presence in the spectrum of a light nucleon made of an ordinary 
quark and a conjugate quark |J. 

This result can be easily understood when we realize that between the onset fi a and the 
saturation point the distribution of the phase of the determinant becomes indistinguish- 
able from a flat one: the averaged determinant is no longer a real and positive quantity as 
it should be. The free energy density defined through the average of the modulus differs 
from the complete one by means of the phase contribution ([]]). With N independent 
configurations the best we can do is to evaluate (e t ^ A )\\ with an error of the order of 
0(^t=). If the contribution of the phase is relevant it has to be exponentially small with 
the system volume but we can not appreciate its contribution unless N = 0(e v ). 

Except for the first, central and last coefficient (constrained to be real and positive from 
eigenvalues symmetries) the same scenario is valid for all the coefficients in the fugacity 
expansion. The coefficients are canonical partition functions (at fixed number density) 
and, in perfect analogy with the grand canonical partition function, the modulus of their 
average practically coincide with the average of their modulus. 

From this considerations we can understand why the Glasgow prescription gives the 
same result as the modulus: unless we have a huge statistics (exponential with the system 
volume) all the possible definitions for a real and positive Z point toward the theory with 
the modulus of the Dirac determinant |J (see for a high statistics study of 2 4 ). The 



MDP algorithm is not plagued from the same problems since one first integrates over the 
gauge fields and then over the matter field. In some sense in this case one avoids the sign 
problem summing over all possible gauge fields exactly. The main conclusion is that the 
moderate optimism that was inspired by the GCPF calculations in the last years has to 
be considered ill-founded. 
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